function [H, RhoExact, EgyExact, mu1] = ham2d5( N, V, mu, h, beta, ...
 isRhoComput, isEgyComput, isMuUsed, NSite );
% Construct the finite difference matrix
ntotal = N^2;

[x,y] = ndgrid( (0:(N-1))/N, (0:(N-1))/N );
nx = N;
ny = N;
display('Generate a 2-d tight binding matrix');
H = sparse(nx*nx, nx*nx);
for i1 = 1 : nx
  for j1 = 1 : ny
    pos = nx*(i1-1) + j1;
    pos1 = nx * (i1-1) + mod(j1,ny)+1;
    pos2 = nx * (i1-1) + mod(j1-2,ny)+1;
    pos3 = nx * (mod(i1,nx)) + j1;
    pos4 = nx * (mod(i1-2,nx)) + j1;
    H(pos, pos1) = -1;
    H(pos, pos2) = -1;
    H(pos, pos3) = -1;
    H(pos, pos4) = -1;
    H(pos, pos) = 4;
  end
end
H = H / h^2;
for i1 = 1 : nx
 for j1 = 1 : ny
   pos = nx*(i1-1) + j1;
   H(pos, pos) = H(pos, pos) + V(j1,i1);
 end
end

if( isRhoComput )
  display('Diagonalizing');
  lc = 8;
  opts.disp = 0;
  % [EF, D] = eig( full(H) );
  [EF, D] = eigs( H, NSite+40, 'sm', opts );
  D = diag(D);
  SD = sort(D)
  if( beta == inf )
    ind = find( D < SD(NSite+1) );
    RhoExact = 2d0 * EF(:,ind) * EF(:,ind)';
    EgyExact = sum(diag( RhoExact * H ));
    RhoExact = reshape(diag(RhoExact), N, N );
    % mu1 = (SD(NSite)+9*1d-9);
    mu1 = (SD(NSite)+SD(NSite+1))/2;
    display('Gap');
    min(abs(mu1-SD))
    pause
    return;
  end
  

  numelectron=@(mux)sum(2./(1+exp(beta*(D-mux))) );
  mu1=fzero(@(mux)numelectron(mux)-2*NSite,mu);
  display('Gap');
  min(abs(mu1-SD))
  % mu1 = SD(NEig)
  pause

  if(isEgyComput)
    display('Constructing density matrix');
    if( isMuUsed == 1 )
      RhoExact = 2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF';
      EgyExact = sum(diag( RhoExact * H ));
    else
      RhoExact = 2d0 * EF * diag( 1./ (exp(beta*(D-mu1))+1d0) ) * EF';
      EgyExact = sum(diag( RhoExact * H ));
      % RhoExact = reshape(diag(RhoExact), N, N );
    end
  else
    display('Constructing density function');
    if( isMuUsed == 1 )
      RhoExact = reshape(2d0 * abs(EF).^2 * ( 1./ (exp(beta*(D-mu))+1d0) ), N, N );
      EgyExact = [];
    else
      RhoExact = reshape(2d0 * abs(EF).^2 * ( 1./ (exp(beta*(D-mu1))+1d0) ), N, N );
      EgyExact = [];
    end
  end
else
  RhoExact = [];
  EgyExact = [];
end



